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Abstract 

A linear method for inverting noisy observations of the Radon transform 
is developed based on decomposition systems (needlets) with rapidly decay- 
ing elements induced by the Radon transform SVD basis. Upper bounds 
of the risk of the estimator are established in (1 < p < oo) norms for 
functions with Besov space smoothness. A practical implementation of the 
method is given and several examples are discussed. 

1 Introduction 

Reconstructing images (functions) from their Radon transforms is a fundamental 
problem in medical imaging and more generally in tomography. The problem is 
to find an accurate and efficient algorithm for approximation of the function to 
be recovered from its Radon projections. In this paper, we consider the problem 
of inverting noisy observations of the Radon transform. As in many other inverse 
problems, there exists a basis which is fully adapted to the problem, in particular, 
the inversion in this basis is very stable; this is the Singular Value Decomposition 
(SVD) basis. The Radon transform SVD basis, however, is not quite suitable for 
decomposition of functions with regularities in other than L^-related spaces. In 
particular, the SVD basis is not quite capable of representing local features of 
images, which are especially important to recover. 

The problem requires a special construction adapted to the sphere and the 
Radon SVD, since usual tensorized wavelets will never reflect the manifold struc- 
ture of the sphere and will necessarily create unwanted artifacts, or will concentrate 
on special features (such as ridgelets...). 

Our idea is to design an estimation method for inverting the Radon trans- 
form which has the advantages of maximum localization of wavelet based methods 
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combined with the stabihty and computabihty of the SVD methods. To this end 
we utihze the construction from [22] (see also [15]) of locahzcd frames based on 
orthogonal polynomials on the ball, which are closely related to the Radon trans- 
form SVD basis. As shown in the simulation section the results obtained are quite 
promising. 

To investigate the properties of this method, wc perform two different stud- 
ies. The first study is of theoretical kind and investigates the possible losses (in 
expectation) of the method in the 'minimax framework'. This principle, fairly 
standard in statistics, consists in analyzing the mathematical properties of esti- 
mation algorithms via optimization of their worst case performances over large 
ensembles of parameters. We carry out this study in a random model which is 
also well known in statistics, the white noise model. This random model is a toy 
model well admitted in statistics since the 80's as an approximation of the 'real' 
model on scattered data. It is proved, for instance in [2] that the regression model 
with uniform design and the white noise model are close in the sense of Le Cam's 
deficiency -which roughly means that any procedure can be transferred from one 
model to the other, with the same order of risk-. This model has the main ad- 
vantage of avoiding unnecessary technicalities. In this context we prove that over 
large classes of functions (described later), our method has optimal rates of con- 
vergence, for all the LP losses. To our knowledge, most of the parallel results are 
generally stated for losses, as in [18] for example, very few (if any) consider 
losses while it is a warrant for instance that the procedure will be able to detect 
small features. Again, the problem of choosing appropriated spaces of regularity 
in this context is a serious question, and it is important to consider the spaces 
which may be the closest to our natural intuition: those which generalize to the 
present case the approximation properties shared by standard Besov and Sobolev 
spaces. We can also prove that our results apply for ordinary Besov spaces. 

In the case p > 4 we exhibit here new minimax rates of convergence, related 
to the ill posedness coefficient of the inverse problem ^ along with edge effects 
induced by the geometry of the ball. These rates are interesting from a statistical 
point of view and have to be compared with similar phenomena occurring in other 
inverse problems involving Jacobi polynomials (e.g. Wicksell problem), see [14]. 

Our second study of the performances of our procedure is performed on sim- 
ulations. Since in practical situations scattered data are generally observed, we 
carried out our simulation study in the scattered data model. We basically com- 
pared our method to the SVD procedure -since it is the most commonly studied 
method in statistics- and the simulation study consistently predicts quite good 
performances of our procedure and a comparison extensively in favor of our al- 
gorithm. One could object that it is a rather common opinion that 'one should 
smooth the SVD'. However, there are many ways to do so (for instance, we men- 



2 



tion a parallel method, employing a similar idea for smoothing out the projection 
operator but without using the needlet construction and in a no-noise framework, 
which has been developed by Yuan Xu and his co-authors in [28, 29, 30]). Ours 
has the advantage of being optimal for at least one point of view since we are able 
to obtain the right rates of convergence in Lp norms. 

The paper is structured as follows. In Section 2 we introduce the model and the 
Radon transform Singular Value Decomposition. In Section 3 we give the class 
of linear estimators built upon the SVD. We also give the needlet construction 
and introduce the needlet estimation algorithm. In Section 4 we establish bounds 
for the risk of this estimate over large classes of regularity spaces. Section 5 is 
devoted to the practical implementation and results of our method. Section 6 is 
an appendix where the proofs of some claims from Section 3 are given. 

2 Radon transform and white noise model 
2.1 Radon transform 

Here we recall the definition and some basic facts about the Radon transform (cf. 
[12], [20], [16]). Denote by B'' the unit baU in R'^, i.e. B"^ = {x = {xi, . . . , xa) e 
R'^ : |a;| < 1} with \x\ = {J2t=i^iV^^ and by S"'"^ the unit sphere in R'^. The 
Lebesgue measure on B'^ will be denoted by dx and the usual surface measure on 
grf-i (Iq-^x) (sometimes we will also deal with the surface measure on S'^ which 
will be denoted by daa). We let \A\ denote the measure \A\ — J^dx ii A G B^ as 
well as \A\ = J^da{x) if A C S"'"^ 

The Radon transform of a function / is defined by 

Rf{e, s)^ j ^^^^ f{s0 + y)dy, 9 e s e [-1, 1], 

where dy is the Lebesgue measure of dimension d — 1 and 9^ = {x G R'^ : {x, 9) — 
0}. With a shght abuse of notation, we will rewrite this integral as 

/ J{y)dy. 

J{y,9)=s 

It is easy to see (cf. e.g. [20]) that the Radon transform is a bounded linear 
operator mapping L^(S'^, dx) into (s*^"^ x [—1, 1], ci/x(^, s)), where 

ds 

dfi{9,s) = da{9)j^—^^y^^. 



3 



2.2 Noisy observation of the Radon transform 

We consider observations of the form 



dY{9, s) = Rf{9, s)dii{9, s) + edW{9, s), 



(2.1) 



where the unknown function / belongs to L^(i?^, dx). The meaning of this equation 
is that for any s) in L^(§'^~^ x [—1, 1], diJi{9, s)) one can observe 



The goal is to recover the unknown function / from the observation of Y . As 
explained in the introduction, this model is a toy model, fairly accepted in statistics 
as an approximation of the 'real model' of scattered data. The study is carried 
out in this setting to avoid unnecessary technicalities. 

Our idea is to devise an estimation scheme which combines the stability and 
computability of SVD decompositions with the superb localization and multiscale 
structure of wavelets. To this end we utilize a frame (essentially following the 
construction from [15]) with elements of nearly exponential localization which is 
compatible with the SVD basis of the Radon transform. This procedure is also to 
be considered as a first step towards a nonhnear procedure especially suitable to 
handle spatial adaptivity since real objects frequently exhibit a variety of shapes 
and spatial inhomogeneity. 

2.3 Polynomials and Singular Value Decomposition of the 
Radon transform 

The SVD of the Radon transform was first established in [5, 6, 17]. In this regard 
we also refer the reader to [20, 28]. In this section wc record some basic facts 
related to the Radon SVD and recall some standard definitions which will be used 
in the sequel. 

2.3.1 Jacobi and Gegenbauer polynomials 

The Radon SVD bases are defined in terms of Jacobi and Gegenbauer polynomials. 
The Jacobi polynomials Pn°''^\ n>0, constitute an orthogonal basis for the space 




^{Rf,(l))^ + eW^. 



Here = J (f){9, s)dW{9, s) is a Gaussian field of zero mean and covariance 




4 



l],WaA't)dt) with weight w^^/sit) = (1 - t)"{l + tf , a,f3> -1. They are 
standardly normahzed by Pi"'^\l) = and then [1, 10, 25] 



where 



.(a,/J) ^ 2"+/^+^ r(n + a + l)r(n + /3 + l) 

(2n + Q; + /3 + l)r(n + l)r(n + Q; + /3 + l)' ^ ' ^ 

The Gegenbauer polynomials are a particular case of Jacobi polynomials, tra- 
ditionally defined by 

^"^'^ = (A + V2)/ "'"''''"'''^^^' ^ > 

where by definition (a)„ = a(a + 1) . . . (a + n — 1) = ^ is readily seen that 

C„^(l) = = and 

(2.3) 



2.3.2 Polynomials on B'^ and § 

We detail the following well known notations which will be used in the sequel. 
Let n„(K'^) be the space of all polynomials in d variables of degree < n. We 
denote by P„(M'^) the space of all homogeneous polynomials of degree n and by 
Vn(M'^) the space of all polynomials of degree n which are orthogonal to lower 
degree polynomials with respect to the Lebesgue measure on B"^. Vo is the set of 
constants. We have the following orthogonal decomposition: 



n„(M'^) = 0v, 



fc=0 



Also, denote by 7in(]R°') the subspace of all harmonic homogeneous polynomials 
of degree n and by 7Y„(§°'~^) the restriction of the polynomials from 7in{^'^) to 
Let nji(S'^~^) be the space of restrictions to S'^"^ of polynomials of degree < n on 
M.^. As is well known 

n 



m=0 



(the orthogonality is with respect of the surface measure da on S" ) 
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Let Yi^i, 1 < i < Nd^i{l), be an orthonormal basis of Hi{S''' i.e. 

Jsd-1 

Then the natural extensions of YJ^j on S'^ are defined by Yi^i{x) — \x\''Yi^i[-^) and 
satisfy 

/ Yi,i{x)Yi~i^dx = [ r-^-^ / YiM)yiA^da{Odr 

J B'^ Jo Jsd-i 

= f r'-^''-' ^ Yi,{OYiA^)da{Odr = 

For more details we refer the reader to [8]. 

The spherical harmonics on S"^"^ and orthogonal polynomials on B"^ are natu- 
rally related to Gegenbauer polynomials. The kernel of the orthogonal projector 
onto Hni^'^'^) can be written as (see e.g. [24]) if Na_i{n) is the dimension of 

Nd-l(:n) 

E YUOYiM- ^J'.^'J.^ C^ii^^^))- (2-4) 

The "ridge" Gegenbauer polynomials Cn'^{{x,^)) are orthogonal to n„_i(i?'^) in 
l^^B*^) and the kernel Ln{x,y) of the orthogonal projector onto V„(-B'') can be 
written in the form (see e.g. [21, 28]) 

Lnix, y) = ^ /' Ci'\{x. 0)C'JMy, 0)d<7{0 (2.5) 

_ {n + l),-i f Ci'\M)Cil\{y,^)) 

The following important identities are valid for "ridge" Gegenbauer polynomi- 
als: 

r uW2) 

/ Cf^((e,^))Cf^((r?,x))rfx = -^C„^/^((e,77)), C,veS''-\ (2.6) 

and, for x e B'^,r) e §'^"\ 

C'JM^,x))C'J\{^,rj))da{0 = \s'-'\C'J\{7j,x)), (2.7) 

1 

see e.g. [21]. By (2.5) and (2.7) 
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/3,1,0 /4,4,1 




/l4,4,l /l4,12,0 

Figure 1: A few radon SVD basis elements (Low quality figure due to arXiv 
constraint) 

and again by (2.5) 

/ Ln{x, OLniy, Od(^iO = (2?^ + d)Ln{x, y). 
2.3.3 The SVD of the Radon transform 

Assume that {Yj^j : 1 < < Nd~i{l)} is an orthonormal basis for 7{i{S'^~^). Then 
it is standard and easy to see that the family of polynomials 

fk,i,i{x) = {2k+dY/^Pf''^''^'~'\2\x\^-l)Yi,{x), 0<l<k, k-l = 2j, l<z< N,_^{l), 
form an orthonormal basis of Vfc(-B'^), see e.g. [8]. On the other hand the collection 

9k,iA(^, s) = [hi'/'Y'^\l-sy'~''>/'Ct^\s)YiM, k>0,l>0,l<t< N,.r{l), 

is obviously an orthonormal basis of L^(§'^^-'^ x [— 1,1], dfi{6,s)). 

Figure 1 displays a few fk,i,i and illustrates their lack of localization. 

The following theorem gives the SVD decomposition of the Radon transform. 

Theorem 2.1. For any f G L^(i?'^) 

A:>0 0<l<k,k-l=0{mod2) l<i<Na-iil) 
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and for any g e l2(S''"^ x [-1, 1], s)) 

R*9 = ^^k {9,9k,i,i)i,fk,i,i- (2.9) 

fc>0 0<l<k,k-l=0{mod2)l<i<Nd-i{l) 

Furthermore, for f e 1?{B'^) 

f = J2^^' E E w^9k,i,i),fk,i,- (2.10) 

fc>0 0<i<fc,A;-Z=0(mod2) l<i<iVd_i(0 

/n ^/le a6ove identities the convergence is in and 

nd—d—l nd—d—1 
\2 _ ^_Jj_ _ ^ ^ i'-'^+l ('9 11'! 

*^ (A; + l)(A; + 2)...(A; + d-l) {k + l)a-i^ ' ^' ' 

Remark : Observe that iik>0, 0<l<k, k — 0(mod2) , and 
l<i<Nd-i{l),thenR*fk,i,i^O. * 

For the proof of this result, we refer the reader to [20] and [28]. We will only 
show in the appendix that has the value given in (2.11). 



3 Linear estimators built upon the SVD 

3.1 The general idea 

In a general noisy inverse model 

dYt = Kfdt + edWt, 

where 7^ is a linear operator mapping / G H i— >■ Kf G K, and IH and K are 
two Hilbert spaces, the SVD yields a family of linear estimators via the following 
classical scheme. 

Suppose K has an SVD 

Kf ^^(Tm{f,em)el,, /GH, 

m 

where {cm} and {e*^} arc orthonormal bases for H and K, respectively, and Kcm — 
arnC^ and K*e*^ = amG-m with K* being the adjoint operator of K. We also assume 
that — ^ 0. Then, if 7^ 0, 

j e*^dYt - j Kf- e*^dt + e J e^^dW^ ^ Jf- K*e*^dt + e J e^^dWt 
= (^m I femdt + 6 e*^dWt 
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and hence 

— / el^dY, = (/, e^) + — / el^dWt. (3.1) 
In going further, suppose that {(pi} is a tight frame for H. Therefore, for any / € H 

f ^^Oii(f)i, ai^{f,(j)i). 
I 

We can represent 0/ in the basis {e^}: 
and hence 

m 

On account of (3.1) this leads to the estimator 

fN = Y.&icl>i with ai = J2^r— fe*mdYt, (3.2) 

where A?" is a parameter. By (3.1) we have 

= E + E / ^IdW^ = ai + Zi, 

J 

m m 

where Zi has a normal distribution -^(0,^^(7/")^^). In this scheme the factors 

-5-, which are inherent to the inverse model, bring instability by inflating the 
variance. 

The selection of the frame {cpi} is critical for the method described above. The 
standard SVD method corresponds to the choice (pi — e;. This SVD method is 
very attractive theoretically and can be shown to be asymptotically optimal in 
many situations (sec Dickcn and Maass [7], Mathc and Pereverzev [19] together 
with their nonlinear counterparts Cavalier and Tsybakov [4], Cavalier et al [3], 
Tsybakov [26], Goldenschluger and Pereverzev [11], Efromovich and Kolchinskii 
[9]). It also has the big advantage of performing a quick and stable inversion of 
the operator K. However, while the SVD bases are fully adapted to describe the 
operator they are usually not quite appropriate for accurate description of the 
solution of the problem with a small number of parameters. Although the SVD 
method is suitable for estimating the unknown function / with an L^-loss, it is 
also rather inappropriate for other losses. It is also restricted to functions which 
are well represented in terms of the Sobolev space associated to the SVD basis. 
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Switching to an arbitrary frame {(pm}, however, may yield additional instability 
through the factors (7™)^'s. 

Our idea is to utilize a frame {4>i} which is compatible with the SVD basis {cm}, 
allowing to keep the variance within reasonable bounds, and has elements with 
superb space localization and smoothness, guaranteeing excellent approximation 
of the unknown function /. In the following we implement the above described 
method to the inversion of the Radon transform. We shall build upon the frames 
constructed in [22] and called "needlets" . 



3.2 Construction of needlets on the ball 

In this part we construct the building blocks of our estimator. We will essentially 
follow the construction from [22]. 

3.2.1 The orthogonal projector Lk on Vfc(-B'^) . 

Let {fk,i,i} be the orthonormal basis of Vfc(-B'') defined in §2.3.3. Denote by the 
index set of this basis, i.e. = {(/, i) : < I < k, I = A;(mod2), < i < Na-i{l)}. 
Also, set u — d/2 — 1. Then the orthogonal projector of L^(S'^) onto Vfe(-B'') can 
be written in the form 

Lkf = j f{y)Lk{x,y)dy with Lk{x,y) = ^ fk,i,i{x)fk,i,i{y). 

i,ieTk 

Using (2.4) Lk{x,y) can be written in the form 

Lk{x,y) (3.3) 

= {2k+d) Yl pr'\2\^\' - iupr^^\2\y\' - En.(^)n.(n) 

l<k,k-l=0i2) i 
' ' l<k,k-l=0{2) 

Another representation of Lk{x,y) has already be given in (2.5). Clearly 

/ Lm{x,z)Lk{z,y)dz = 6m,kLm{x,y) (3.4) 

and for / e l\B'^) 

f^Y^^kf and 11/11^ = E 11^^/112 = E(W,/>- (3.5) 

fc>0 k k 



fI \v\ 
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3.2.2 Smoothing 

Let a e C°°[0, oo) be a cut-off function such that < a < 1, a{t) = 1 for f e [0, 1/2] 
and supp a C [0, 1]. Wc next use this function to introduce a sequence of operators 
on h'^{B'^). For j > write 



^jfi^) = y^(^(^^Lkf{x) = I Aj{x,y)f{y)dy with Aj(a;, y) = V a^^^^L 
Also, we define Bjf = Aj+i/ — Ajf. Then setting b{t) = a{t/2) — a{t) we have 
Bjfix) = J2b{^)Lkf{x) = j^B,{x,y)f{y)dy with B,{x,y) = ^^2^)^^ 



X 



k 



Evidently, for / e 1?{B- 

k 

k 

and 



{Ajf,f)^J2^{^)(^'^f^f)^\\f\\l (3-6) 



hm \\A,f-f\\2= lim ||(Ao + y]i?^)/-/||2 = 0. (3.7) 



J— >0O J— >0O 

m=0 



An important result from [22] (see also [15]) asserts that the kernels Aj{x,y), 
Bj{x,y) have nearly exponential localization, namely, for any M > there exists 
a constant cm > such that 

\AJx,y)\,\BJx,y)\ < Cm , =, x,y e B'^, (3.8) 

I A ,y)u\ A ,y)\- ^^(i + 2^d(x,y))^yw^yw^' ' ^ ^ 

where Wj{x) — + a/I — \x\'^, |a;p = |a;|^ = a;^, and 

d 

d{x, y) = Arccos ((x, y) + a/1 - - ?/) = ^iVi- (3-9) 



i=l 



The left part (1) of Figure 2 illustrates this concentration: it displays the 
influence of a point x to the value of Bjf at a second point y^, namely the values 
of Bj{x,yo) for a fixed i/q and j — 4. This influence peaks at yo and vanishes 
exponentially fast to as soon as one goes away from Uq. The central part (2) of 
Figure 2 shows the modification of the concentration when j is set to a large value 
(j = 6). The right part (3) of Figure 2 shows the lack of concentration of Bj when 
the cut-off function a used is far from being C°°. The resulting kernel still peaks 
at yo but the value of Bjf at yo is strongly influenced by values far away from yo. 
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j 4 6 4 

a C°° C°° non smooth 

(1) (2) (3) 

Figure 2: Smoothing kernel Bj{x,yo) for a fixed yo for (1) a C°° cut-off function 
a with j = 4 , (2) for the same a with j = Q and (3) for a non smooth cut-off 
function a with j = 4 (Low quahty figure due to arXiv constraint) 

Remark : At this point it is important to notice the following correspondence 
which will be used in the sequel. For S'l = {{x, z) G x ]R+, + z'^ = 1}, we 
have the natural bijection 

X & B"^ X = {x, v^l — and d{x,y) = dgd_{x,y), 

where d^d is the geodesic distance on §;jj^. 



3.2.3 Approximation 

Here we discuss the approximation properties of the operators {^j}- We will show 
that in a sense they are operators of "near best" polynomial L^-approximation. 
Denote by En{f,p) the best L^-approximation of / G IJ'{B'^) from n„, i.e. 

E4f,p) = mf \\f-P\\,. (3.10) 
Estimate (3.8) yields (cf. ([22, Proposition 4.5]) 

/ \Aj{x,y)\dy<c*, x G 5^ j > 0, 

where c* is a constant depending only on d. Therefore, the operators Aj are (uni- 
formly) bounded on L^{B'^) and L°^(i?'^), and hence, by interpolation, on L^(i?'^), 
1 < p < oo, i.e. 

\\AjW<c*\\f\\„ fei7{B'). (3.11) 

On the other hand, since a(t) = 1 on [0, 1/2] we have AjP = P for P G n2i-i. 
We use this and (3.11) to obtain, for / G IJ'{B'^) and an arbitrary polynomial 



12 



P e n2.-i, 

= \\f-P+P-A,f\\, < ||/-P||^+p,(P-/)||, < (l+c*)||/-P||, = K\\f-P 

Consequently, ||/ — Aj/||p < KE2]-i{f,p). In the opposite direction, evidently, 
Ajf e Ha. and hence ^2^(/,p) < ||/ - Ajf\\p. Therefore, for / e L^'(P"'), l<p< 
oo, 

^2.(/,p) < 11/ - AjfWp < KE^-.{f,p). (3.12) 

These estimates do not tell the whole truth about the approximation power of Aj. 
It is rather obvious that because of the superb localization of the kernel Aj{x,y) 
the operator Aj provides far better rates of approximation than E2j-i{f, oo) away 
from the singularities of /. 

In contrast, the kernel Sj{x.y) = ^o<fc<2J ^k{x,y) of the orthogonal projector 
Sj onto II23 is poorly localized and hence Sj is useless for approximation in L^, 
Py^ 2. This partially explains the fact that the traditional SVD estimators perform 
poorly in Z^-norms when p ^ 2. 

3.2.4 Splitting procedure 

Let us define 



Cjix,y) = ^Ja(^^^Lmix,z)) and Dj{x,y) = ^ Jb(^^^Lm{x,z). 

m m 

Note that Cj and Dj have the same localization as the localization of Aj, Bj in 
(3.8) (cf. [22]). Using (3.4), we get the desired splitting 

Aj{x,y)= Cj{x,z)Cj{z,y)dz (3.13) 

and 

Bj{x,y)^ f Dj{x,z)Dj{z,y)dz. (3.14) 

J Ed 

Obviously z ^ Cj{x, z)Cj{z, y) is a polynomial of degree < 2-'+^ and z ^ P'j{x, z)Dj{z, y) 
is a polynomial of degree < 2^^^. The next step is to discretize the kernels Aj{x, y) 
and Bj{x, y). 

3.2.5 Cubature formula and discretization 

To construct the needlets on B'^ we need one more ingredient - a cubature formula 
on B'^ exact for polynomials of a given degree. 
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Recall first the bijection between the ball B'^ (equipped with the usual Lebesgue 
measure) and the unit upper hemisphere in M'^+^: 

S^ = {(x,y), XGM^ 0<y<l, \x\l + y^^l} 

equipped with da the usual surface measure. 

T : (x, y) e 1-^ X e M'^ 

and 

T-^ : X eR'^ ^ X ^ {x, ^l- e 
Applying the substitution T one has (see e.g. [27]) 

/ F{x,y)da{x,y)^ [ F{x, y/l - \x\^) (3.15) 
and hence for / : M'' i— > M 



L 



dx 



f{x)ydad{x,y) = f{x)Jl- \x\l = / f{x)dx. (3.16) 

{x,y)eSf JxeBd " VI ~ Fid "'S'* 

Therefore, given a cubature formula on S"' one can easily derive a cubature formula 
on B'^. Indeed, suppose we have a cubature formula on §^ exact for all polynomials 
of degree n + 1, i.e, there exist Xn C §^ and coefficients 0;^ > 0, ^ e Xn, such that 



/ 

J si 



P{u)da{u) = Yl ^1^(0 VP e n„+i(M'^+i). 



If P e n„(M'^) then P{x)y e n„+i(K'^+^) and hence 

Vc^|P(e)yr^- / P{x)yda^ [ P{x)dx. 

iexu + 

Thus the projection Xn of Xn on is the set of nodes and the associated coefficients 
given by = a/1 — ■C^cj^ induce a cubature formula on B^ exact for n„(M'^). 
The following proposition follows from results in [22] and [27]. 

Proposition 3.1. Let {B{C,i,p) : i & 1} be a m,aximal family of disjoint spherical 
caps of radius p = t2^^ with centers on the hemisphere SlJ.. Then for sufficiently 
small < T < 1 the set of points Xj = {d '■ i ^ 1} obtained by projecting the 
set {^i : i e /} on B'^ is a set of nodes of a cubature formula which is exact for 
n2j+2 {B'^) . Moreover, the coefficients cu^^ of this cubature formula are positive and 
~ W^j(6)2"^''^. Also, the cardinality 4xj ~ 
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3.2.6 Needlets 

Going back to identities (3.13) and (3.14) and applying the cubature formula 
described in Proposition 3.1, we get 

Aj{x,y)^ Cj{x,z)Cj{z,y)dz =y2uj^Cj{x,^)Cj{y,^) and 
{x,y) = / Dj{x,z)Dj{z,y)dz = ^uj^Dj{x,i)Dj{y,^). 



We define the father needlets 0^,^ and the mother needlets V'j,^ by 

(pj,i{x) = ^Cj{x,^) and i/jj,^ix) = ^Dj{x,^), ^ e Xj, J > 0- 
We also set ■0-i,o — jbS] X-i — {0}- Prom above it follows that 

Therefore, 

Ajfix) = / y)f{y)dy = = Yl "j.^'^j.^ "j.^ = (/' '^j.^)- 

(3.17) 

and 

= / Bj{x, y)f{y)dy = if^ ^u)^J,i = J2 /^i-^^J.^ ^i^^ = (/' V'i.^)- 

(3.18) 

By (3.17) and (3.6) we have 
and hence 



Il0,,dl2 < 1- (3.19) 

From (3.5) and the fact that ^j>ob{t2~^) = 1 for t e [l,oo), it readily follows 
that 
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and taking inner product with / this leads to 

ii/ii2 = EEK/'^^-.^)i'' 

which in turn shows that the family {ij^j,^} is a tight frame for Ij^{B'^) and conse- 
quently 

Il^.,dl2 > Il^.,dl2, i-e. Il^,,dl2<l- (3.20) 

Observe that using the properties of the cubature formula from Proposition 3.1, 
estimate (3.8) leads to the localization estimate (cf. [22]): 

2jd/2 

10,- .(x) I, < Cm , VM > 0. (3.21) 

Nontrivial lower bounds for the norms of the needlets are obtained in [15]. More 
precisely, in [15] it is shown that for < p < oo 

Ifellp ~ ll'/'.Jp ~ , i^Xr (3.22) 

We next record some properties of needlets which will be needed later on. For 
convenience we will denote in the following by hj^^ either 0^^^ or 

Theorem 3.2. Let 1 < p < oo and j > 1. The following inequalities hold 

J] ||%^||^ < c2^'(*/2+(p/2-2)+) (3 23) 

J] < cj2^*/' if p = 4, (3.24) 

and for any collection of complex numbers {d^}^^^. 

II E d^h^^P <c{Yl imh^di)'^'. (3.25) 
^ex^- ^ex^- 

Here c> is a constant depending only on d, p, and r. 

To make our presentation more fluid we relegate the proof of this theorem to 
the appendix. 
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3.3 Linear needlet estimator 



Our motivation for introducing the estimator described below is the excellent ap- 
proximation power of the operators Aj defined in §3.2.2 and its compatibility with 
the Radon SVD. We begin with the following representation of the unknown func- 
tion / 

/ ~ ^ fk,l,i) fk,l,ii 

where the sum is over the index set {{k, : k > 0,0 < I < k,l = A;(mod2), 1 < 
i < Nd^i{l)}. Combining this with the definition of Aj we get 

where 

= ^J2ikj,i{fJk,i,i) = J2iki,iY^W'>^9k,i,X ^J^ikhy^ J 9k,i,iR{fW- 

Here 7^'|. = {fk,i,i, (pjM) can be precomputed. 

It seems natural to us to define an estimator fj of the unknown function / by 

fj^Yl (3-26) 
«.,« = E^&^ / 9k,i,idY. (3.27) 

k,l,i ^ 

Here the summation is over {{k, l,i) : < k < 2\0 <l < k,l = /c(mod 2), 1 < i < 
A^d_i(Z)} and j is a parameter. 

Some clarification is needed here. The father and mother needlets, introduced 
in §3.2.6, are closely related but play different roles. Both (pj^^ and ipj^^ have su- 
perb localization, however, the mother needlets have multilevel structure 
and, therefore, are an excellent tool for nonlinear n-term approximation of func- 
tions on the ball, whereas the father needlets are perfectly well suited for linear 
approximation. So, there should be no surprise that we use the father needlets for 
our linear estimator. 

Furthermore, even if the needlets are central in the analysis of the estimator, 
the estimator fj can be defined without them. Indeed, 



where 



fo = Y^H^khy 1 9k,i,idY <f>j,^ 



^exj k,l,i 
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as all the sum are finite, their order can be interchanged, yielding 



fj = HYkj ^^'^'^^ E ^kh^i.^ = E ^ / 9k,i,idY Ajfk,i,i 

k,l,i CSXj k,l,i 

and thus the estimator is obtained by a simple componentwise multiplication on 
the SVD coefiicients 

fj = E / gk,i,idY fk,i,i. 

However, as will be shown in the sequel, the precise choice of this smoothing 
allows to consider hp losses and precisely because of the locahzation properties 
of the atoms this approach will be extended using a thresholding procedure in a 
further work, using this time the mother wavelet. 

4 The risk of the needlet estimator 

In this section we estimate the risk of the needlet estimator introduced above in 
terms of the Besov smoothness of the unknown function. 

4.1 Besov spaces 

We introduce the Besov spaces of positive smoothness on the ball as spaces of 
L*'- approximation from algebraic polynomials. As in §3.2.3 we will denote by 
£■„(/, p) the best L*'- approximation of / e 'U'{B'^) from n„. We will mainly use 
the notations from [15]. 

Definition 4.1. [15] Let 0<s<oo, l<p<oo, and < q < oo. The space 5^'° 
on the ball is defined as the space of all functions f e L^(S'^) such that 

\^\bP° = (j2'^n'E^{f,p)y-y^'' < oo ifq<oo, 

77,/ 

n>l 

and I/Ib^.o = sup^^i n^En{f,p) < oo if q — oo. The norm on B^'^ is defined by 

ll/llB-» = ll/llp + l/li.;:0- 

Remark : From the monotonicity of {En{f,p)} it readily follows that 

ll/llB|;0~ll/llp+(E(2''^2.(/,l>)r)'^' 
j>0 
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with the obvious modification when q — oo. ★ 

There are several different equivalent norms on the Besov space B^'^. 

Theorem 4.2. With indexes s,p,q as in the above definition the following norms 
are equivalent to the Besov norm 



« = 11(2^11/ -A,/y,>oii^., 
(ii) M2{f)^\\f\\r,+mn\Bjfu>i\\i., 

Proof. The equivalence A/i(/) ~ 11/ Us"." is immediate from (3.f2). 

To prove that A/"2(/) ~ Afiif), we recall that Bj = (Aj+i-Aj) (see §3.2.2) and 
hence \\Bjf\\p < which readily implies A/'2(/) < aA/'i(/). 

In the other direction, we have 

oo oo 
l=j 1=3 

Assuming that A2(/) < oo we have = ai2~''^ with {ai} e l'^. Hence 

oo oo oo 

1=3 l=j l=j 

and by the convolution inequality {I3j} G l'^. Therefore, A/i(/) < aA/2(/). 

For the equivalence Nz{f) ~ II/IIb^.O) see [15, Theorem 5.4]. □ 

4.1.1 Comparison with the "standard" Besov spaces 

The classical Besov space B^^^{B'^) is defined through the L^-norm of the finite 
differences: 

^hf{x) = {f{x + h)- /(x))l^eB<il^+^gBci 

and in general 

N / \ 

^hfi^) = ^.^B'^^.+NkeB^J^^-lf-'rAfix + kh). 

k=0 ^ ^ 

Then the A^th modulus of smoothness in is defined by 

a;,^(/,t) = sup||A^/||„ t > 0. 

\h\<t 
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For 0<s<A^, l<p<oo, and < g < oo, the classical Besov space B^^ is 
defined by the norm 

= 11/11. + - (f;r<(/,2-^-)]')'^' 



j=Q 

with the usual modification for o = oo. It is well known that the definition of B?. „ 
does not depend on N as long as s < N [13]. Moreover, the embedding 

b;,, q (4.1) 

is immediate from the estimate En{f,p) < cWp (/, 1/n) [13]. 

4.2 Upper bound for the risk of the needlet estimator 

Theorem 4.3. Let l<p<oo,0<s<oo, and assume that f e -Bp ^ ^^^^ 
11/11^.0^ < M. Let 

be the needlet estimator introduced in §5.5, where J is selected depending on the 

param,eters as described below. 

1. // M2-''(*+'^) ~ e whenp^ oo, then 

E||/ - /jIIoo < CooM^£^VlogM/£. 

2. If M2--^' ~ £2-^(d-2/p) y^^g^ 4 < p < oo, t/ien 

^ (d-2/p)p sp 

- < CpM-+d-2/P£«+d-2/p, 

where when p — 4 there is an additional factor ln{M/e) on the right. 

3. IfM2-^' - £2'^('^-V2) when l<p<A, then 



^ (d-l/2)p ap 



Remarks : 



• It will be shown in a forthcoming paper that the following rates of conver- 
gence are, in fact, minimax, i.e. there exist positive constants ci and C2 such 
that 

~ sp sp 

sup _ inf E||/ - /ll^ > ci max{£^+rf-2/p, £«+d-i/2}^ 

||/|Ls,o <M /estimator 

sup _ inf E||/-/||oo > C2e^y/logl/e. 

ll/ILs,o <M/ estimator 
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• The case p = 2 above corresponds to the standard SVD method which in- 
volves Sobolev spaces. In this setting, minimax rates have aheady been 
cstabhshcd (cf. [7], [19] [4], [3], [26], [11], [9]); these rates are e^+^ . Also, 
it has been shown that the SVD algorithms yield minimax rates. These re- 
sults extend (using straightforward comparisons of norms) to LP losses for 
p < 4, but still considering the Sobolev ball {||/||oo,o < M} rather than 

2,oo 

the Besov ball {||/||^s,o < M}. Therefore, our results can be viewed as an 



extension of the above results, allowing a much wider variety of regularity 
spaces. 

• The Besov spaces involved in our bounds are in a sense well adapted to our 
method. However, the embedding results from Section 4.1.1 shows that the 
bounds from Theorem 4.3 hold in terms of the standard Besov spaces as well. 
This means that in using the Besov spaces described above, our results are 
but stronger. 

• In the case p > 4 we exhibit here new minimax rates of convergence, related 
to the ill posedness coefficient of the inverse problem ^ along with edge 
effects induced by the geometry of the ball. These rates have to be compared 
with similar phenomena occurring in other inverse problems involving Jacobi 
polynomials (e.g. Wicksell problem), see [14]. 




4.3 Proof of Theorem 4.3 




(4.2) 



Now from 



dY = RfdiJi + edW 



we have 





/ fk,i,i dx + e Zk,i,i 



f R*9k,i,i dx + e Zk,i,i 



and hence 
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On account of (3.27) this leads to 

JBd .,, 



Here the summation is over {{k,l,i) : < k < 2\0 < I < k,l = k{mod 2),1 < i < 
Nci-i{l)}. Since Z^^i^i are independent iV(0, 1) random variables, Z^^g ~ A^(0, a|^) 
with 

4. = E hiifj^ < ^ S (4.3) 

with c = ((i/27r)'^~^. Here we used that {fkn} is an orthonormal basis for and 

hence E.,M 170/ = ll0.,dli<l- 

Prom (3.26) fj — X^^g^^. ctj,^0j,^and using (4.2) we have, whenever 1 < p < oo, 

- ^11^ < 2^-^11/ - A,f\\l + E\\A,f - f^} 

< 2^-\cM^2-^^^ + n\Ajf - IX) (4.4) 

and, for p = oo, 

E||/ - ^-lloo < 11/ - ^.-/lloo + n^if - ^ Hoc 

<cM2-^-^ + E||^,/-^-|U. (4.5) 
On the other hand, using inequality (3.25) of Theorem 3.2 we obtain, if 1 < p < oo, 

and hence 

Eii^.-/ -f3\\i<cY. nzi^^rUidi < c{e2^^'-'^/'r UjM. (4-6) 

where we used that B\Zj^^\P < c{e2^^'^-^^/'^)P. Similarly, for p = oo, 

W^jf - fjWoo = II 5^(tti,e - ^m)<Pm\\oo < cmaxjaj-^ - Sj,^|||0j,dloo 

and hence 

n^jf - fjWoo < cE{max|Z,-^|||0,-^||oo} 

< 082^^"-'^^ max \\(i)j,^\\^^/2h^ < 082^"^. (4.7) 
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For the second inequality above we used Pisier's lemma: If Zj ~ A'"(0, of)-, Oj < cr, 
then 

E( sup \Zj\) < a^/2hg^. 

l<j<N 

We also used that max^^^^. H^j^^Hoo < c2^(^+'^y'^^ which follows by inequality (3.23) 
of Theorem 3.2. 

Combining (4.5) and (4.7) we obtain, for p — oo, 

E||/-^||oo<c{M2-^-^ + £2^'^V7} 

and if M2-^(^+'^) ~ s, then 

E||/ - ^ lloo < cM^£^ VlogM/£. 

Similarly, combining estimate (3.23) of Theorem 3.2 with (4.6) and inserting the 
resulting estimate in (4.4) we obtain in the case 4 < p < oo 

= c{MP2-^'"P + £P2^'(*-2)}. 
If M2-^' ~ £2j(d-2/p) i^his yields 

^ (d-2/p)p sp 

E||/ - fjWl < CM s + d-2/pirs + d-2/p_ 

Accordingly, for p = 4 we combine inequality (3.24) with (4.6) and insert the result 
in (4.4) to obtain 

- ^11^ < c{MP2-^'P + (£2^(<^-i)/2)P_^-2i*/2| 
= c{MP2-^'P + j{e2^^'^-^/y} 
and if M2-i' ~ e2^id-im this yields 

Ell/ 11^ < cM ^+d-2/P e log M/e. 
Finally, if 1 < p < 4 as above we obtain using (3.23), (4.6), and (4.4) 
E||/ - fjWl < c{MP2-^'P + (^s2^(d-i)/2y2:idp/2y 
= c{MP2-^'P+{e2^^'^-^/y}. 
So, if M2-^' ~ £2^(d-i/p)^ then 

^ (d-l/2)p sp 

E||/ - fjWl < cM-+'i-l/2£«+d-l/2. 

This completes the proof of the theorem. □ 
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Figure 3: Simplified CAT device 



5 Application to the Fan Beam Tomography 
5.1 Radon and 2d Fan Beam Tomography 

We have implemented this scheme for d = 2 in the radiological setting of Cormack[5] 
This case corresponds to the fan beam Radon transform used in Computed Axial 
Tomography (CAT). As shown if Figure 3, an object is positioned in the middle of 
the device. X rays are sent from a pointwise source S{9i) located on the boundary 
and making an angle 9i with the horizontal. They go through the object and are 
received on the other side on uniformly sampled array of receptors R{9i, 92). The 
log decay of the energy from the source to a receptor is proportional to the integral 
of the density / of the object along the ray and thus one finally measures 



for 9 = 9i ~ 92 and s = sin ^2- The device is then rotated to a different angle 
9i and the process is repeated. Note that d9-^j^^ is nothing but the measure 
corresponding to the uniform d9id92 by the change of variable that maps (6*1, ^2) 
into {9, s). 

The Fan Beam Radon SVD basis of the disk is tensorial in polar coordinates: 



h,i,i{r,9) = {2k + 2y/'P^'''\2\r\' -l)\r\%,,{9), < I < k, k-l = 2j, l<t<2, 





[-1,1], 



where P°'' is the corresponding Jacobi polynomial, and Yii{9) = cicos{l9) and 
'^1,2(9) = ci sm{9) with cq = and q = ^ otherwise. The basis of x [—1, 1] 
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has a similar tensorial structure as it is given by 



guUO.t) = [hk]-^'\l-t'f'^Cl{f)Yu{e), A;>0, />0, l<z<2, 

where Cl is the Gegenbauer of parameter 1 and degree k. The corresponding 
eigenvalues are 

VkTl 

In this paper, we have considered the theoretical framework of the white noise 
model. In this model, we assume that we have access to the noisy "scalar product" 
J gk,i,idY, that is to the scalar product of Rf with the SVD basis gki,i up to a i.i.d. 
centered Gaussian perturbation of known variance e^. This white noise model is 
a convenient statistical framework closely related to a more classical regression 
problem with a uniform design on 9i and 62, which is closer to the implementation 
in real devices. In this regression design, one observe 



Yi,,i, = Rf \^27i - — j ,sin27r— j +ei„,„ h < Ni, i2 < N2 

where A''i and A^2 gives the discretization level of the angles 9i and 6*2 and ej^^jj is 
an i.i.d. centered Gaussian sequence of known variance a^. Note that this points 
are not cubature points for the SVD coefficients. The correspondence between the 
two model is obtain by replacing the noisy scalar product J gk,i,idY with by the 
corresponding Riemann sum 

^ ^ il=0 12=0 \ \ 1 ^/ ^/ 

and using the calibration = l{N\ x A^2)- It is proved, for instance in [2] that 
the regression model with uniform design and the white noise model are close in 
the sense of Le Cam's deficiency -which roughly means that any procedure can 
be transferred from one model to the other, with the same order of risk-. The 
estimator jj defined in the white noise model by 



h — y 9k,i,i(iy jk,i,i — — — j gk,i,iay jk,i 

is thus replaced in the regression model by 



- ^ a (^) v"^ ^ iii) y/k + 1 

fj = 2^ — \Rf, 9k,i,i) fk,i,i = 2^ 7r~/= 9k,i,i) fk,i 
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5.2 Numerical results 



To illustrate the advantages of the hnear needlet estimator over the linear SVD 
estimator, we have compared their performances on a synthetic example, the clas- 
sical Logan Shepp phantom[23], for different norm and different noise level, and 
for both the white noise model and the regression model. The Logan Shepp phan- 
tom is a synthetic image used as a benchmark in the tomography community. It is 
a simple toy model for human body structures simplified as a piecewise constant 
function with discontinuities along ellipsoids (see Figure 6). This example is not 
regular in a classical sense. Indeed, it belongs to but not to any 5^'° with 
s>l. 

To conduct the experiments, we have adopted the following scheme. Denote 
by / of the Logan Shepp function presented above, its decomposition in the SVD 
basis fk^i^i up to degree k — 512 has been approximated with an initial numerical 
quadrature x valid for polynomial of degree 4 x A; = 2048, 

(/, fk,l,i) - X] ^{rifii)fi'^i^ di)fk,l,i{ri, Oi) = Ck,l,i 
{ri,0i)ex 

and used this value to approximate the original SVD coefficients of R{f), the 
noiseless Radon transform of /, 

{R{f),fk,i,i) ^ Kck,i,i- 

In the white noise setting, for all k < ko = 256, a noisy observation J gk,i,idY 
is generated by 

9k,LidY ~ XkCkjA + eWk,i,, 



where e is the noise level and Wk.i.i a lid sequence of standard Gaussian random 
variables. Our linear needlet estimator fj of level J — \og2{k^) is then computed 

as ^ 

k<ko,l,i ^ 

while the linear SVD estimator f^s of degree k^ is defined as 

fks^y2 + T-Wk,i,i)fk,i,i- 

Ak 

We also consider the naive inversion up to degree ko which is equal to /jf^ . The 
estimation error is measured by reusing the initial quadrature formula, 

11/ -/r^ E ^inA)\f{r^A)-f{r^A)Y. 
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In the regression setting, we have computed the values of the Radon transform 
Rf of / on a equispaced grid for the angles 6i and 62 specified by its sizes A^^i and 
N2 using its SVD decomposition up to A; = A; = 512. We have then defined the 
noisy observation as 

with ej^^jj an i.i.d. centered Gaussian sequence of known variance cr^. The esti- 
mated SVD coefficients are obtained through the Riemann sums 

ii=0 12=0 ^ ^ ' ' 

We plug then these values instead of the J gk,i,idY in the previous estimators. 

For each noise level and each norm, the best level J and the best degree K 
has been selected as the one minimizing the average error over 50 realizations of 
the noise. Figure 4 displays, in a logarithmic scale, the estimation errors ||/ — 
in the white noise model plotted against the logarithm of the noise level e. It 
shows that, except for the very low noise case, both the linear SVD estimator and 
the linear Needlet estimators reduce the error over a naive inversion linear SVD 
estimate up to the maximal available degree /cq. They also show that the Needlet 
estimator outperforms the SVD estimator in a large majority of cases from the 
norm point of view and almost always from the visual point of view as shown in 
Figure 6. The localization of the needlet also 'localizes' the errors and thus the 
"simple" smooth regions are much better restored with the needlet estimate than 
with the SVD because the errors are essentially concentrated along the edges for 
the needlet. Remark that the results obtained for the regression model in Figure 5 
are similar. We have plotted, in a logarithmic scale, the estimation errors against 
the logarithm of the equivalent of the noise in the regression a^/ (Ni x A^2) with 
Ni = N2 = 64 and various a^. Observe that the curves are similar as long as o"^ is 
not too small, i.e. as long as the error due to the noise dominate the error due to 
the discretization. As can be seen both analysis do agree. This confirms the fact 
that the white noise model analysis is relevant for the corresponding fixed design. 

A fine tuning for the choice of the maximum degree is very important to obtain 
a good estimator. In our proposed scheme, and in the Theorem, this parameter 
is set by the user according to some expected properties of the unknown function 
or using some oracle. Nevertheless, an adaptive estimator, which does not require 
this input, can already be obtained from this family, for example, by using some 
aggregation technique. A different way to obtain an adaptive estimator based on 
thresholding is under investigation by some of the authors. 
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Figure 6: Visual comparison for the original Logan Shepp phantom with 
(Low quality figure due to arXiv constraint) 
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6 Appendix 

6.1 Proof of identity (2.11) 

From [20, p. 99] with some adjustment of notation, we have 

where < I < k and I = k (mod 2). As will be seen shortly Aj^ is independent of I. 

We will only consider the case d > 2 (the case d = 2 is simpler, see [20, p. 
99]). To compute the above integral we will use the well known identity (cf. [25, 
(4.7.29)]) 

Summing up these identities (with indices n,n — 2, . . . ) and taking into account 
that C^{t) = 1, C^{t) = 2X{t), we get 

c'^^\t) = y: I cu.it). (6.1) 

3=0 

This with A = {d — 2)/2 and the orthogonality of the polynomials Cn~'^^^'^{t), 
n > 0, yield 

' ci^\t)ci'-'^/\t){i - t'f-'y^dt = ^ [ct'^^'imi - t^-'^/'dt 

1 A J_i 

l + ^AX) (/ + A)2i-2^7r r(/ + 2A) 



A ' Ar(A)2 (/ + A)r(/ + i)' 

We use this and that C^(l) = §2.3.1) and iS-^-^l = r{{d-iyl) ^° obtain 

27r<^-i k\r{d)l\r{d-2) 2'^-'^7r r{l + d-2) 



>^k r(^)r(^)r(^ + d)r(z + d-2)(^_2)r(^)' r(z + i) 

2'-'^n'^ r{d)T{d-2) 1 



r(¥)r(^) (d - 2)T{^Y {k + i)d-i ■ 

The doubling formula for Gamma-function says: r(2^) = ^^^^r{z)r{z + 1) (see 
e.g. [25]) and hence 

r{d)r{d-2) = {d-i){d-2)r{d-2r = ——{d-i){d-2)r[-^) . 

We insert this in (6.2) and then a httle algebra shows that A^ = (fc+'i)^^ • ^ 
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6.2 Proof of Theorem 3.2 

For the proof of estimates (3.23)-(3.24) we first note that by (3.3) 

1 



and we need an upper bound for 2~^''y',^, ^ To this end, 

we will use the natural bijection between B"^ and considered in the remark in 
§3.2.2. Thus for x e B'^ write x = {x, ^1 - \x\^) e S^. Let p = (0, 1) be 
the "north pole" of S'^. For ^ e Xj we denote by Ssd(^,p) is the geodesic ball 

on S'^ of radius p centered at i.e. B^d{^,p) := {£ G S*^ : d^d{x^p) < p}, where 
dsd{x,p) = Arccos (^1 — |a;p) = Arccos {x,p) is the geodesic distance between 
X, p. Using that - |^|| < \{u,p) - {i,p)\ < dsd{i,u) < p for « e and 
p = t2^^ < 2"^ (see Proposition 3.1), it follows that 

< . . . = .^ yueB^4i,p). 



On the other hand, we have \B^d{t p)\ = IS'^'^] J^{sm9Y-^d9 > p'^|S"'-i|| 

with |S'^~-'^| = r{d/2) • above and the fact that the balls {5sd(|, p)}gg^ 

are disjoint to obtain 

< 2-'^ Y ^ f . \^ da{u) 



1 , ... , fsin^)-^-! 



- (2-^ + (w,P))'- ^ ^ ~ ' ' Jo (2-^ + cos^)' 

Jo (2 J + cos^)'^ Jo (2 ^ + ty J2-j 

This yields estimates (3.23)-(3.24). 

We now turn to the proof of estimate (3.25). We will employ the maximal 
operator M.i (t > 0), defined by 

Mtf{x) := sup [ \f{y)\'dy) \ x e B^, (6.3) 

B3X V|-D| Jb / 

where the sup is over all balls B C B'^ with respect to the distance d{-,-) from (3.9) 
containing x. It is easy to show that (see §2.3 in [15]) the Lebesgue measure on B'^ is 
a doubling measure with respect to the distance d{-,-). Hence the general theory 
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of maximal operators applies. In particular, the Fefferman-Stein vector- valued 
maximal inequality is valid: IfO<p<oo,0<g<oo, and < t < min{p, q} then 
for any sequence of functions {fiy}u on B"^ 

oo oo 

IKE \MtM-rf% < cii(E (6.4) 

Denote by B{^,r) the projection of i?sd(^,r) onto B'^, i.e. B{^,r) := {x & B'^ : 
d{x,^) < r}. By [15, Lemma 2.5], we have 

{MtlBi^,r)){x) > c(l + e e 5', < r < TT. (6.5) 

It is easy to see (cf. [15]) that 

|S(e,p)|~2-^'^(2-^ + V^r^)~2-^-''iy,(0, UXr (6.6) 

Also, we let !£■ := i^Ie denote the L^-normalized characteristic function of E G 
B'^. Then (3.21) and (6.6) imply 

ll^j,dlp~ C&Xj- (6.7) 

Now, pick < t < 1 and M > ((i+ l)/t. From (3.21) and (6.5) it follows that 

K^{x)\ < c{MtlBii,,)){x), X e B". (6.8) 

Using this, the maximal inequality (6.4), and (6.7) we obtain 

II E ^^^idp < c|| E Mt{d^lB{^,p))\\p < c\\ E d^^B{^,p)\\p 



This completes the proof of (3.25). □ 
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